function TradeLevel
%% LOAD DATASET
load SST_panel_3; 
% Set data size 
n = 136; m = n; s = 1;
nobs = n; n = nobs; m = n; dim = 5;
% Define dependent variable
Y =  TRADE(s:s-1+n,s:s-1+n)';
% Define independent variables
X = zeros(n,m,dim);
X(:,:,1 ) = LNDIST(s:s-1+n,s:s-1+n)';
X(:,:,2 ) = BORDER(s:s-1+n,s:s-1+n)';
X(:,:,3 ) = COMLNG(s:s-1+n,s:s-1+n)';
X(:,:,4 ) = COLONY(s:s-1+n,s:s-1+n)';
%X(:,:,5 ) =   OPEN(s:s-1+n,s:s-1+n)';
%X(:,:,5 ) = COMFRT(s:s-1+n,s:s-1+n);

%X(:,:,5 ) =   OPEN_WTO(s:s-1+n,s:s-1+n)';
X(:,:,5 ) = COMFRT_WTO(s:s-1+n,s:s-1+n)';

         
 XX = cell(dim,1); for k = 1:dim, XX{k} = X(:,:,k); end; X = XX;

% Pseudo Poisson
pmle = [-.75; .37; .38; .08; .38];
% GMM
[gmm, ~, se_gmm, ~, ~, cond, it] = Newton(@GMM2,pmle,Y,X,{eye(dim)},1),











function [likelihood score Hessian mVar Se asyvar] = GMM2(psi,Y,X,W,pow)
% dimensions
[n m] = size(Y  );  % sample size
[K D] = size(psi); % number of regressors, number of GMM estimators (different sets of moments) 
J = zeros(D,1); for d = 1:D, J(d) = size(W{d},1); end % number of moments J for each D
nn = nchoosek(n,2); mm = nchoosek(n-2,2); rho = nn*mm;
% variable definitions
I = cell(D,1); for d = 1:D, I{d} = zeros(n,m); for k = 1:K, I{d}   = I{d} + X{k}*psi(k,d); end; end
T = cell(D,1); for d = 1:D, T{d} = zeros(n,m); for k = 1:K, T{d}   =            exp(I{d}); end; end
V = cell(D,K); for d = 1:D,                    for k = 1:K, V{d,k} =           T{d}.*X{k}; end; end
% weight matrix inversion
for d = 1:D, W{d} = inv(W{d}); end
% construction of empirical moments
q  = zeros(D    ,1); 
p  = zeros(K    ,1);
%p1 = zeros(pow  ,1); % continuous
%p2 = zeros(1    ,1); % binary
%p3 = zeros(pow-1,1); % interactions
Q  = zeros(K    ,1);
S  =  cell(D    ,1); for d = 1:D, S{d} =zeros(J(d),1); end
H  =  cell(D    ,1); for d = 1:D, H{d} =zeros(J(d),K); end
% variance pre-allocation
mVar    = cell(D,1); for d = 1:D, mVar{d}    = zeros(J(d)   ,J(d)); end
Upsilon = cell(D,1); for d = 1:D, Upsilon{d} = zeros(n,m,J(d),n,m); end
Z = X{1}.* X{2};
for i1 = 1:n,
    for j1 = 1:m,
        if j1==i1; continue; end
        for i2 = i1+1:n,
            if i2==j1; continue; end
            for j2 = j1+1:m,
                if j2==i1 || j2==i2, continue; end
                % SCORE VECTOR
                % generalized residual
                for d = 1:D
                    q(d) = Y(i1,j1)*Y(i2,j2)*T{d}(i2,j1)*T{d}(i1,j2)-T{d}(i1,j1)*T{d}(i2,j2)*Y(i2,j1)*Y(i1,j2);
                end
                % instrumental variables for each estimator 
                for k = 1:K,
                    p(k,:) = (X{k}(i1,j1)-X{k}(i1,j2))-(X{k}(i2,j1)-X{k}(i2,j2));
                end
                % moment conditions
                for d = 1:D, 
                    S{d} = S{d} + p*q(d)/rho; 
                end
                
                % MOMENT VARIANCE
                for d = 1:D,
                    if d==1, 
                        Upsilon{d}(i1,j1,:,i2,j2) = p*q(d);
                        Upsilon{d}(i2,j1,:,i1,j2) = p*q(d);
                        Upsilon{d}(i1,j2,:,i2,j1) = p*q(d);
                        Upsilon{d}(i2,j2,:,i1,j1) = p*q(d);  
                    end
                end
                
                % JACOBIAN MATRIX
                % derivative of generalized residual
                for d = 1:D,
                    for k= 1:K,
                        Q(k,d) = Y(i1,j1)*Y(i2,j2)*T{d}(i2,j1)*T{d}(i1,j2)*(X{k}(i2,j1)+X{k}(i1,j2))-T{d}(i1,j1)*T{d}(i2,j2)*Y(i2,j1)*Y(i1,j2)*(X{k}(i1,j1)+X{k}(i2,j2));
                    end
                end
                % Jacobian matrix
                for d = 1:D, 
                    if d==1, 
                        H{d} = H{d} + p*Q(:,d)'/rho; 
                    end
                end
                
                
            end
        end
        
    end
end

likelihood = cell(D,1); for d = 1:D, likelihood{d} = zeros(1,1); end
score      = cell(D,1); for d = 1:D, score{d}      = zeros(K,1); end
Hessian    = cell(D,1); for d = 1:D, Hessian{d}    = zeros(K,K); end

for d = 1:D, 
    likelihood{d} = -  S{d}'*W{d}*S{d}; 
    score{d}      = -2*H{d}'*W{d}*S{d};
    Hessian{d}    = -2*H{d}'*W{d}*H{d};
    
    xi = 4*mean(mean(Upsilon{d},5),4)*n/(n-2)*m/(m-3); 
    for k1 = 1:J(d),
        for k2 = 1:J(d),
            mVar{d}(k1,k2) = mean(mean(xi(:,:,k1).*xi(:,:,k2)))*n/(n-1);   
        end
    end
    
    Avar = inv(H{d}'*W{d}*H{d})*(H{d}'*W{d}*mVar{d}*W{d}*H{d})*inv(H{d}'*W{d}*H{d});
    Se(:,d) = sqrt(diag(Avar)/(n*m)); asyvar = Avar/(n*m);
end

likelihood_temp = likelihood{1}; likelihood = likelihood_temp;
score_temp      = score{1};      score= score_temp;
Hessian_temp    = Hessian{1};    Hessian = Hessian_temp;
%Se_templ        = Se{1};         Se = Se_temp;
mVar_temp       = mVar{1};       mVar = mVar_temp;
